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Abstract 
Q 

r \ I Recent studies on the phase-space dynamics of a one-dimensional Lennard-Jones fluid reveal 

^ ■ the existence of regular collective perturbations associated with the smallest positive Lyapunov 

exponents of the system, called hydrodynamic Lyapunov modes, which previously could only be 
identified in hard-core fluids. In this work we present a systematic study of the Lyapunov exponents 
QQ ' and Lyapunov vectors, i.e. perturbations along each direction of phase space, of a three-dimensional 

'sj- ■ Lennard-Jones fluid. By performing the Fourier transform of the spatial density of the coordinate 



part of the Lyapunov vector components and then time-averaging this result we find convincing 



o 

QQ I signatures of longitudinal modes, with inconclusive evidence of transverse modes for all studied 



densities. Furthermore, the longitudinal modes can be more clearly identified for the higher density 



^ , values. Thus, according to our results, the mixing of modes induced both by the dynamics and the 

Cu I dimensionality induce a hitherto unknown type of order in the tangent space of the model herein 

studied at high density values. 
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I. INTRODUCTION 

In recent years it has been possible to study in detail the underlying chaotic dynamics 
of many-particle systems, thus finding interesting connections between their microscopic 
dynamics and the observed macroscopic behavior. For example, in the case of static prop- 
erties, the largest Lyapunov exponent (LLE), which measures the dynamical instability of 
phase-space trajectories to infinitesimal perturbations in the initial conditions, is related to 
the Kosterlitz-Thoules transition temperature in a system of coupled rotators with nearest- 
neighbor interactions [l|], as well as to first [2] and second-order phase transitions |3| in 
particle systems with long-range interactions. For dynamical properties, the sudden change 
in the gradient of the LLE against energy, which corresponds to the transition from weak to 
strong chaos, can be detected in the macroscopic behavior of a Brownian particle coupled 
to the system ^. It has also been related to transport coefficients of fiuid systems with 
continuous potentials in non-equilibrium situations |5|. In equilibrium, a relation has also 
been proposed p, although for this case arguments both against |7| and in favor [8] have 
been advanced. Nevertheless, these last results entail interesting possibilities, since they sug- 
gest that dynamical instability is at the origin of macroscopic transport phenomena. Two 
different methodologies have related the Lyapunov spectrum (LS), defined as the sorted set 
of Lyapunov exponents (LEs) which give the exponential rate of expansion or contraction of 
nearby trajectories along each independent component of the phase space, to the transport 
coefficients of simple fiuids J9|, [lOj. However, they involve setting up non equilibrium simula- 
tions or locating special phase space trajectories. Furthermore, neither of these approaches 
considers the perturbations in phase space underlying the LEs, which have attracted a lot 
of attention in recent years. Therefore it would be desirable first to attain a more thorough 
understanding of the phase-space perturbations associated to the LEs of the most general 
type of dynamical models employed to study fiuid systems before attempting the construc- 
tion of a general theory that could relate, with enough confidence, the dynamical instability 
with the macroscopic behavior as described by the transport coefficients. 

For some particular many-particle systems, e.g. hard-sphere fiuids, the theory of LEs 

n 

is highly developed [ll!]. Nevertheless, their most interesting features have been discovered 
by means of molecular dynamics simulations, which revealed that the slowly-growing and 
decaying perturbations associated with the non-vanishing LEs closest to zero are related to. 



and in some cases may even be represented as, almost exact periodic vector fields coherently 
spread out over the physical space with well defined wave vectors, together with an almost 
precise stepwise structure of the considered LEs |12| . From their similarity to the classical 
modes of fiuctuating hydrodynamics, these perturbations have been named Lyapunov modes 
(LMs). Since their discovery [l3|, much work has been done to understand their origin and 
possible relevance. Some analytical approaches that have been advanced to understand 



these small LEs include random-matrix dynamics [ij, llSJ] , periodic orbit models 16| , and 



Q[ 






kinetic theory |17l . [iSj. All these approaches have met with only a very limited success, 
since none of them has been able to satisfactorily describe the mode dynamics in all the 
known situations accessible to simulations. A more serious drawback is that neither has 
been capable to extract the transport coefficients from the LMs, or to relate them, in a 
simple way, to hydrodynamic fiuctuations. 

Until recently the existence of these LMs could only be verified in one, two and three- 

n 

dimensional hard-sphere fluids pjj; their existence in the case of atomic fluids with soft 
interactions (both attractive and repulsive) remained controversial [19|, since certain fea- 
tures, such as the step structure in the LS, dissappear in soft-potential systems. Hence, to 
extend the concept of LMs beyond the realm of hard-core systems, a more generic deflnition 
has to be adopted, which involves both the appearance of a sharp peak of low wave number 
in the spatial Fourier spectrum of LVs corresponding to the LEs closest to zero and of a 
minimum in the average spectral entropy, deflned in Sec. IIVBI of the present work, in that 
same LE region. Employing this more generic deflnition convincing evidence of the exis- 
tence of the so deflned LMs was obtained for some soft-potential systems, namely in fluids 
with a one- dimensional (ID) Lennard- Jones (LJ) [20|, where the aforementioned deflnition 
of LMs and the spectral analysis techniques needed to detect them were developed, and 
two-dimensional (2D) Weeks- Chandler- Anderson (WCA) interaction potentials J2l|. More 
recently, the existence of LMs has been also verifled in lattices of coupled Hamiltonian and 



dissipative maps 22|, |23|, l2J] , as well as in Fermi-Pasta-Ulam anharmonic oscillator lattices 
251] . Thus all information so far available seems to indicate that the Hamiltonian structure, 
conservation laws, and translational invariance are not necessary conditions for the existence 
of LMs. Nevertheless, since there is no theoretical scheme that can predict the existence 
of LMs, the discovery and characterization of LMs has to be done in a case-by-case basis. 
For example, the introduction of damping in a system of coupled Hamiltonian maps does 



not destroy the LMs, whereas the addition of that same damping to a system of coupled 
circle maps wipes off the LMs [23)]. Therefore, although LMs are present in the aforemen- 
tioned systems, their existence is in no way guaranteed in other situations not studied so far, 
such as in three-dimensional Hamiltonian fluid systems with both attractive and repulsive 
interactions. 

In this paper we perform a study of the dynamical instabilities of a three-dimensional 
(3D) atomic fluid interacting with the full LJ potential under the simulation conditions most 



frequently encountered in molecular dynamics studies [26|. We assess the validity of various 



dynamical indicators previously proposed in the literature for our particular system. 



The 
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main result of our work is that, although the spectral analysis methods proposed in Ref. 
can be readily applied to our system, the dimensionality greatly enhances the mixing among 
perturbations that was already present in the case of the ID LJ gas, thus rendering both 
the detection and characterization of the collective perturbations even more problematic. 

This paper is organized as follows: In Sec. [TTl we present a survey of the model, as well 
as the computational details needed to obtain our results in the phase space of our system 
as well as a short account of the relevant theory necessary to study the perturbations of 
the phase space. Sec. [Till describes the results for the complete LS, its dependence on the 
particle density, and provides evidence of the mixing among perturbations. In Sec. IIVI 
we show the results of the spectral analysis methods applied to the proposed dynamical 
indicators. Section |V] is devoted to discuss and analyze the results reported in the previous 
sections. In the Appendix we report some preliminary results for the temporal correlations 
of the phase space perturbations. We present our conclusions in Sec. |VT1 

II. THE MODEL AND SIMULATION DETAILS 
A. Phase-space dynamics 



The L J potential f/^^ (r) between particles i and j is given by 



t^..W=4e 



(J\i2 /cr\6 
r / Vr / 



(1) 



where r = rij is the distance between particles i and j; here a and e are the LJ atom 
diameter and the strength of the interparticle interaction, respectively. The actual interac- 
tion potential employed in our simulations is the spherically truncated and shifted (STS) 



potential, which can be written as 

y r > Tc, 

where Vc = 2.5a is the cut-off radius at which the potential is truncated in order to save 
computer time. It is important to notice that the force derived from the above potential, as 
well as its derivative, are not continuous at the truncation point r^. Therefore the following 



potential, first proposed in Ref. [27|, can be employed 



r > Tc, 

with C2 = 4e[6((T/rJi2 _ ^[a/r^f] and f/^ = 4e[-7((T/r,)i2 ^ 4((x/rj6]_ T^is potential has 
a continuous derivative in the truncation point r^. Nevertheless, for the ID LJ model it 
has been shown, by adopting the above expression for the potential, that the qualitative 
)ehavior of the LMs is not greatly affected by the discontinuity of the STS potential at r^ 
20(]. In the next section it will be shown explicitly that this same phenomenology holds for 
our 3D LJ fluid. The complete Hamiltonian of our system is then 

N 2 

1=1 ' i<j 

with {pi} being the momenta of the atoms and {mj} their corresponding masses. In all 
our simulations we took the masses of the atoms equal to unity, i.e. rrii = m = 1\/ i, with 
(7 = 1, and e = 1 as well. The Hamiltonian is then written in terms of the reduced variables 
r* = r/a, p* = pi/ (me) , and m* = mi/m. 

The initial configuration for all simulations was set up from a fee lattice on a square cell 
of sides L* = L* = L* = L*,. The initial momenta were drawn from a Maxwell-Boltzmann 
distribution. Then the 6A^ equations of motion were numerically integrated by means of 
the Verlet leap-frog algorithm with periodic boundary conditions and the minimum image 
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convention applied in all directions. A time-step of At* = At (e/m) /a = 0.001 was used 
in all simulations, with an equilibration period of 100 000 steps at constant temperature 
obtained by a uniform rescaling of the velocities at each time-step. After the equilibration 
period the system was allowed to evolve at constant total energy for a period of 10^ time- 
steps. The relative energy drift for this number of time-steps is 10~'^-10~^, which is an 



acceptable compromise between accuracy and speed, since there is no systematic drift and 
thus the energy fluctuations are stable for the chosen time-step. All the reported results were 
obtained for systems of A^ = 108 at a supercritical reduced temperature of T* = fc^T/e = 1.5, 
where k^ is the Boltzmann constant. The reduced densities p* = pa^ were taken within 
the range p* G [0.01,0.5]. The critical temperature and density for the LJ fluid with STS 
potential are T* = 1.085 and p* = 0.317 [2&\; thus the employed values for these variables 
ensure that the system state is far away from the two-phase region in the p* vs T* phase 
diagram, being a homogeneous fluid. Finally, since in the rest of the work reduced variables 
will be exclusively employed, from now on we will drop the asterisk from all symbols without 
risk of confusion. 

B. Tangent-space dynamics 

The phase space trajectory is represented by the variable T{t) = 
(Fi (t) , r2 (t) , . . . , Tm (t)), where Fj (t) = (fj (t) , pj (t)). To study the local dynam- 
ical stability of our system we introduce the Lyapunov vector (LV) as ST^^'^t) = 
(5fS°) (t) , 5r(") (t) , . . . , 6T^;:^ (t)) , where a = 1, . . . , 6Ar and ^fJ") (t) = {6rf (t) , 6pf (t)) 
representing the ith particle contribution to the infinitesimal perturbations of the trajectory 
F (t) along all possible directions (position and momentum axes) of the phase space, thus 
defining the so called tangent space. 



According to Oseledec's multiplicative ergodic theorem [29|] the remote past limit sym- 
metric operator $b(^) = limio^-oo[-^(^;^o) ■-^"'"(^j^o)]^''^^^*"*"^' exists for almost every initial 
condition F (to), where A4(t,to) is the fundamental matrix governing the time evolution of 



the perturbations ST (t) in tangent space as ST{t) = Ai{t,to) ■ 6T{to) [lO|. The set of in- 
stantaneous LEs is defined as \^°'\t) = In A^"\ where A*^"") are the eigenvalues of $b(^)- The 
herein employed standard procedure to compute the LEs consists in periodically reorthonor- 



malizing a set of offset vectors that are time evolved by means of the matrix Ai {t, tg) 30|, l31 |. 
The time-averaged values of the logarithms of the renormalization factors, i.e. (In A'^"))^ are 
the LEs JA*^"-*} and the set of offset vectors right after the reorthonormalization are the 
eigenvectors of $b(^)j which are called backward LVs. The equivalence of the LEs computed 
by means of these two methods can be proved rigorously [32] , but the relation between the 
Oseledec eigenvectors and the LVs obtained via the standard method is more subtle. It is 



known that the backward LVs converge at an exponential rate to the Oseledec eigenvectors 



for the inverse-time dynamics of the original system 33|, |3J] . The latter are obtained as the 
eigenvectors of the far future limit operator $f(t) = \iintg^oo[-M.^{to,t) ■ X(to,^)]^^^^*-*°~*''' 
and are called forward LVs. It is to be noted that recently it has been possible to obtain from 
the intersection of the embedded subspaces spanned by the eigenvectors of $b(^) (backward 
LVs) and $f(t) (forward LVs) the so called characteristic LVs, which are independent of the 
norm and do not form an orthogonal basis [35]. A somewhat similar algorithm has also been 



proposed in Ref. [36||. The aforementioned schemes are based on ideas already discussed 
long ago 371], but have only been recently proposed because their implementation is by no 
means a simple task from a computational point of view, which is why they only have been 
tested in simple ID systems. Thus it is not at all clear that their implementation could 
be feasible in the near future for the case of our system. Furthermore, since the purpose 
of the present work is to investigate the possible existence of LMs in a 3D LJ fluid, it is 
important to recall that the discovery of LMs in all studied systems so far has been made 
employing the backward LVs. So it seems that the use of the CLVs is not essential for the 
detection of LMs. Finally, meaningful results are still obtained by means of the backward 
LVs obtained from the standard procedure, such as in the recent characterization of LMs in a 



diatomic system [38| and the discovery of LMs in the XY rotator model [39[; in both studies 
the backward LVs were employed. Therefore, from now on we will refer to the numerically 
computed vectors (backward LVs) as the LVs without confusion. 

In the Hamiltonian case (which we treat here), and also in some special homogeneous non- 



equilibrium situations 40(], the LEs and the corresponding LVs have a symmetry property 
which makes unnecessary to calculate the whole spectrum. In these cases the LS thus 
computed is symmetrical around zero, which means that each LE has a partner that is 



exactly its negative. This is the so called conjugate-pairing rule [41|- Therefore, only 3A^ 
linearized equations for the LVs were simultaneously integrated, along with the 6A^ nonlinear 
equations for the reference trajectory r(t). The initial LVs ^F^^^ (to) in all our simulations 
consisted in a set of 3A^ orthogonal vectors with randomly selected components. 




FIG. 1: LS JA^"^} as a function of the Lyapunov index a computed for a system with the STS 
potential (circles) and the potential defined by Eq. ([3]) (diamonds). The thermodynamic state in 
both cases is defined by p = 0.1 and T = 1.5, with N = 108. 

III. LYAPUNOV SPECTRUM 



As mentioned in the previous section, both the force derived from the STS potential, 
as well as its derivative, are not continuous at Vc- In order to assess the relevance of this 
particular feature in the computation of the properties in the tangent space we perform 
additional simulations employing the potential given by Eq. ([3]). In Fig. [1] we present 
the LS computed from this last potential, as well as that computed by means of the STS 
potential, for T = 1.5, p = 0.1, and A^ = 108. As can be readily appreciated there are no 
significant differences, a result which supports the use of the simpler STS potential. The 
most relevant feature of this figure is that, as in the case of the ID LJ gas, in the smallest 
positive LE region there is no evidence of the stepwise structure that signals the appearance 
of the LMs in the case of hard-core systems. Thus it is plausible that the same mechanism 
that accounts for the absence of the stepwise structure in the LS of the ID LJ gas is also at 
work in the present model. 

In Fig. [2] we present the probability density function for the instantaneous LEs A'^^-'and 
;\(320) ^Qj^ A^ = 108 and p = 0.01. It can be observed that the fluctuations around their mean 
values increase as the Lyapunov index a increases, i.e. when going from a high to a low 
value of the LE. It is to be observed that the fluctuations in the A^^^^-* value are so great that 
the tails in the probability density function overlap with those corresponding to A*^^-*. These 
large fluctuations are certainly a reason why no stepwise structure can be found in the low-a 
region of the LS. Other dynamical indicators that we will present later on are also affected 
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FIG. 2: Probability density function of the instantaneous LEs A'^^ = 0.51 (circles) and A^^^*^' = 
0.003 (squares) for p = 0.01. Filled symbols correspond to the results obtained with the potential 
given by Eq. ([3]). Average taken over 2 x 10^ integration time-steps. Same A^ and T values as in 
Fig. [TJ Vertical dot-dashed lines indicate the values for A'-^) and A*^^^'^^ 



by this behavior of the LS at large a values. In the same figure results for the potential 
given by Eq. (JSj) are also presented; these are quite similar to those obtained with the STS 
potential. Thus it is confirmed that the aforementioned fluctuations are not a result of the 
discontinuity of the STS potential at Vc- Therefore this potential will be used in the rest of 
this work. Other dynamical indicators that we will present later on are also affected by this 
behavior of the LS at large a values. 

In Fig. [3] we show the normalized LS for N = 108 and T = 1.5 computed for the 
employed range of reduced densities. The values of the scalar length of on edge of the cubic 
simulation cell, as well as the values of the LLE corresponding to each value of the reduced 
density, are reported in Table [11 As can be readily appreciated in Fig. [3], the LS {A*^"^} 
corresponding to the lowest density p = 0.01 can be separated into two regions. In both 
of them the LS is a decreasing function of the Lyapunov index a, but in the former the 
decrease is more pronounced than in the latter. This bending of the LS has been observed 
in a quasi ID hard-disk gas [42] as well as in the ID LJ gas 20|, and has been related to the 
separation of two time scales. To properly explain this point it is important to remember 
that each LE indicates a time scale given by its inverse, so the LS can be considered as a 
spectrum of time-scales. The smallest positive LE region of the spectrum is dominated by 
macroscopic time and length scale behavior. On the other hand, the opposite region of the 
LS is dominated by short time scale behavior, such as local collision events. As the density 




FIG. 3: Normalized LS JA^^'/A^^'} as a function of the Lyapunov index a for all employed values 
of the reduced density p and a reduced temperature of T = 1.5, with A^ = 108. 

TABLE I: Length L of the cubic simulation cell along with the values of the LLE for each reduced 
density value. 



p 


L 


Ai 


0.01 


22.104 


1.014 


0.025 


16.287 


1.531 


0.05 


12.927 


2.033 


0.1 


10.260 


2.606 


0.25 


7.560 


3.496 


0.5 


6.000 


4.349 



increases, the collisions increase both in frequency and in the number of particles involved; 
the correlation among them increases, and becomes increasingly difficult to distinguish them 
individually. Thus it is no longer possible to make a separation of time scales, and so the 
LS does not present the aforementioned bending for p > 0.1. Although this explanation 
is plausible for the LS of our 3D LJ fluid, we also point out that the bending depicted in 
Fig. [3]for p = 0.01 (the lowest reduced density value employed) is less pronounced than the 
corresponding lowest density instance of the LS of the ID LJ gas [20]. We interpret this 
fact as a signature that, even for the lowest density case, the effect of the dimensionality 
in the microscopic events is to enhance the correlation among them, and thus reduce the 
separation of time scales that is more evident in the ID case. This explanation will be 
further supported in the next sections. 
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IV. SPATIAL STRUCTURE OF TANGENT SPACE PERTURBATIONS 

A. Spatial Lyapunov vector density 

As discussed in Sec. IIIH the strong fluctuations of the smaller LEs make the perturbations 
in tangent space extremely unstable, rendering any coherent structure that may exist in 
configuration space diflicult to detect. To investigate the possibility that these structures 
(LMs) exist in the case of the 3D STS LJ potential we have to establish a measure of the 
contribution of a given LV ST^°'^ (t) at each point r of the conflguration space, regardless 
of which particle makes the contribution to the magnitude of the chosen LV. We have 
to remember that the ith. particle contribution to the inflnitesimal perturbation ^r*^") (t) 
consists of spatial and momentum components. Since the LMs can be in general considered 
as Goldstone modes resulting from translational invariance in coordinate space [18], we will 
consider only the spatial part of the full perturbation component ST^ (t). Thus, in analogy 
with the deflnition of microscopic density fluctuations 43| , we deflne the spatial LV density 
as 

N 

u(")(r,t) = J]5r;")5(r-r,), (5) 



which was flrst introduced in Ref. 4J] and is the function to be studied afterwards. In 
Fig. Hlwe present a snapshot of a single component u" iz,t) = X]j=i ^z^^ (i" — ^i) of the 
aforementioned spatial density along the z axis of our system. Since the latter is isotropic, 
the results are completely analogous to those corresponding to the Ux and Uy components 
projected along their respective coordinate axes. We can readily appreciate that the spatial 
density corresponding to the LLE is more localized than that corresponding to the LE with 
a = 320, which is the lowest index value corresponding to a LV not related to the space 
and time translational invariance symmetries of the system and the associated conserved 
quantities, the total energy and the total momentum [17|. A point to be noticed is that the 
localization of the LLE is not so clearly deflned and strong as in the ID LJ system [20(| or 



in hard-disk systems [42| . This is another indication that the dimensionality of the system 
has indeed a strong influence on the tangent-space dynamics. 

Next, we proceed to consider the spatial Fourier transform of u*^"-* (r, t), which can be 
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FIG. 4: Snapshot of the spatial LV density component u^ {z, t) along the z axis of the simulation 
box for a LJ fluid of A^ = 108 atoms at a reduced density of p = 0.01 corresponding to the LVs (a) 
a = 1 and (b) a = 320. 



written as 



5(")(k,t) = /"u(°)(r,t)exp(-ik-r)rfr 



N 



5^5rr)exp[-ik-r,(t)] 



i=l 



To proceed further we invoke the static LV density correlation function of Ref. 4^, defined 



as 



C(") (k, t) = G(°) (k, t) G(") (-k, t) 



(6) 



which in the case of our 3D system is a second rank tensor. For our isotropic fluid the 
Cartesian components C/^" (k, t) of C^"^ (k, t) can be written in terms of longitudinal C^"^ 
and transverse C^"'^ static correlation functions as 



C^^J {k,t) = k,KC["^ {k,t) 



flU '"/l' 



^^Uc(")(k,t) 



(7) 



with kfj, = (k/k) . From this last expression it is immediate to obtain the explicit form of 
the longitudinal and transverse static correlation functions as 



Cf ) (k, t) 



^tiu V^i ^) f^ukfj, 



^(cl:HKt)-clzHKt)Kk) 
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To simplify the analysis we will consider a coordinate system such that the wave vector 
k is parallel to the z axis. Then 

AT 

a(°) (k, t) = Y, 5r^) exp \-\Kzi {t)\ . (8) 



i=l 



For homogeneous systems with translational invariance the Fourier transformed components 
of a quantity such as the spatial LV density < Uu k where u = x,y, z, are uncorrelated J45|. 
Thus, our simplification in no way destroys essential information. The spatial Fourier spec- 



,(") 



trumPr(A;„t) = u)r' {h,t) 



,(") 



r(") 



corresponding to each component of the spatial density of 



^Fj can be readily computed by an algorithm for unequally-spaced data points [46|] which 



has been previously applied to the ID LJ '20| and the 2D WCA systems [2l|]. Further- 
more, we observe that the diagonal components of the static LV density correlation function 
correspond to the spatial Fourier spectrum, i.e. ClZ {kz,t) = Pu^ {kz,t). Thus the static 
longitudinal and transverse correlation functions can be obtained from the aforementioned 
power spectra as 






2 

Finally, since averaging over several spatially equivalent directions will improve the statistics, 
we take successively the k vector along the remaining coordinate axes x and y to obtain 
two more sets of longitudinal and transverse correlation functions and then average over all 
sets. The result are longitudinal C^°'\k,t) and transverse C^°'\k,t) correlation functions 
independent of the employed coordinate system. 

At variance with the hard-core systems in which the patterns resembling transverse modes 



do not survive time averaging 19j, the Fourier spectral techniques so far presented have been 
quite successful in the case of the ID LJ system, and so there is a reasonable possibility of 
success in the case of our 3D LJ fluid. 

B. Time instability of instantaneous quantities 

Due to their mutual interaction, the LMs in all soft potential systems are only of finite 
life-time. To investigate their time stability we present in Fig. \5\ the time evolution of two 
quantities associated with the longitudinal static correlation function C^\k,t) for the LV 
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a = 320: the peak wave number /cmax, which indicates the position of the highest peak in the 
spatial Fourier spectrum Pzz{k,t), and the spectral entropy H^"\t) |47|, which measures 
the k distribution properties of the aforementioned spectrum. This last quantity was first 



employed in Ref. 20|] in the study of the ID LJ gas and is defined as 

if (") (t) = - ^ C(") {k, t) In Cj") {k, t) . (9) 

k 

For the ID LJ system these quantities show an intermittent behavior, i.e. large intervals 
of nearly constant low values {off state) axe interrupted by short periods of bursts {on state) 



where they experience large values 20|]. As can be appreciated in Fig. [5](a) for the time 
evolution of the peak wave number fcmax, this behavior is somewhat different in our case; 
although there is an alternation between the on and off states, the mixing among them 
is so great that there are no long-lived time intervals for either state, in sharp contrast 
with the ID case. The instant value of the spectral entropy, presented in Fig. Mjo), has a 
similar behavior as that of fcmax in the sense that it is not easy to identify a correspondence 
between its temporal evolution and that of the on and off states. A virtually identical 
behavior is obtained for fcmax and H^\t)^ but now defined in terms of CJ"^(A;, t) (not shown), 
and for all values of the reduced density. This intermittency in the time evolution of the 
spatial Fourier spectrum of LVs is a typical feature of soft-potential systems. However, 
this behavior is greatly increased in comparison to the ID case. There are two possible 
causes for this difference with the ID results: first, the reduced temperature at which we are 
working is much higher (an order of magnitude) than that employed in Ref. [20|; second, the 
dimensionality of our system is also higher, thus the mixing is easier by this enlargement in 
the phase space, as already explained. Now, in order to average out temporal fluctuations, 
and thus extract useful information about the collective modes, from now on we will study 
the properties of the average spectra (C^^\k,t))^ and (C^\k^t))^^ where {■■■)^ means 
temporal average. Further insight into the temporal dynamics of the LMs can be obtained 



by studying the dynamic LV density correlation function first introduced in Ref. [4^ for 
the ID LJ gas. In the Appendix some preliminary results of the aforementioned dynamical 
correlation function applied to our 3D LJ fluid will be presented. 
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FIG. 5: (a) Time evolution of the peak wave-number fcmax and (b) spectral entropy H^' (t) for 
a Lyapunov index value of a = 320, for A^ = 108, p = 0.01, and T = 1.5. The behavior of these 
quantities corresponding to the transverse correlation function C}^'{k,t) (not shown) is almost the 
same as that depicted in the above figure. 

C. Time-averaged power spectrum 

Figures[6](a) and[6]^b) display, for p = 0.01, (C^"^(A;, t))^ and (C^"^(A;, t))^ respectively. In 
both cases we can readily appreciate that the contribution of the high wave number com- 
ponents is not small. The instantaneous power spectra are not dominated by a single peak; 
rather, several small peaks, which are related to intermediate length scales, are present, 
which in turn make significant contributions to the overall shape of the time averaged cor- 
relation functions for large k values. This feature can be explained by the less pronounced 
time-scale separation than in the ID case, as mentioned in Sec. IIIII in relation to the LS. 
Nevertheless, the highest value of the time averaged correlation functions is always domi- 
nated by certain low-wave-number components; for the longitudinal correlation function we 
observe that the sharp-valued peak in the spectrum corresponds to diminishing /cmax values 
as the Lyapunov index a — > 3A^, i.e. as we go from the region of high LEs to the region of 
low LEs. At this point it is important to notice that this correspondence is not monotonic, 
since the A^max value of the highest peak is attained at a = 200, and the height of the corre- 
sponding peak slightly diminishes, although remains well defined, for higher a values. For 
the transverse correlation function depicted in Fig. Mjo) the results are similar, except for a 
very important feature: the highest value of the spectra is again attained at a = 200, a value 
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FIG. 6: (a) Time averaged longitudinal correlation function (^C^'{k,t)Y and (b) time averaged 
transverse correlation function {C^'{k,t)\ for various values of the Lyapunov index a, with pa- 
rameters N = 108, T = 1.5, and p = 0.01. In the inset the static structure factor S{k) of the 
system is plotted. 



far away of the region corresponding to the lowest LEs, but then vanishes as a ^ 3A^. This 
feature is not consistent with the existence of well-defined transverse LMs, since it would be 
expected that the highest peak in the transverse correlation function should correspond to 
small values of /cmax in the a ~ 3A^ regime. In the inset the static structure factor S{k) 43 1 
for the corresponding thermodynamic state is plotted. The highest peak of this function 
is located at k/2Tx ^ 0.08. Now, the position of the the peak in (^C^\k,t)Y is located at 
^max ~ 0.04. This value suggest that there is no obvious direct connection between the short 
range order of the atoms and the lowest wavevector peak of the longitudinal LV correlation 
function. 

The longitudinal and transverse correlation functions corresponding to the reduced den- 
sity p = 0.5 are presented in Figs. [Tl^a) and[7](b), respectively. A difference with respect to 
the low density results is that the highest peak in each power spectra is broader than in the 
corresponding low-density case. This can be interpreted as an indication that the influence 
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FIG. 7: (a) Time averaged longitudinal and (b) transverse correlation functions for the same 
parameters as in Fig. ([6]), but with a reduced density of p = 0.5. 

of higher wave numbers is stronger than in the low-density case. We further observe that 
the highest peak is located at a higher wave number value than the corresponding peaks 
in Figs. El^a) and (b). Now, although at this density there is a higher degree of spatial 
order in the atoms of the system, the lowest wavevector peak of S{k), which is plotted in 
the inset, has no relation whatsoever with the fcmax value of the time averaged transverse 
LV correlation function. Next we notice that the highest peak of the transverse correlation 
function is reached at a value of the Lyapunov index of a = 200, with a monotonic decrease 
for higher a values. This same result was obtained for the transverse correlation function 
in the low-density regime, although in this case the peak remains well defined, in contrast 
with the result displayed in Fig. [6](b), which is almost a flat spectrum. 

So far our results are consistent with the existence of longitudinal LMs. To quantify the 
properties of these modes in Fig. [8](a) we plot /cmax versus the Lyapunov index a. It is 
clear from this figure that, as a — ;► 3A^, fcmax diminishes. A feature that stands out is that, 
despite the time averaging over a large number of data points (10^ time steps), the obtained 
values remain somewhat noisy. Nevertheless it is clear from the figure that the decrease 
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in the fcmax value as a increases is approximately monotonic, in clear contrast to the ID 
case |[20] in which a sudden change of fcmax from a finite value to zero was interpreted as a 
signature of the separation of time scales mentioned in Sec. IIIII Next, in Fig. [Hl^b) the height 
<^C|^")(/Cjnax)i)). of the highest peak in the time averaged longitudinal correlation function is 
also plotted as function of the Lyapunov index a. On average we can observe a monotonic 
increase of the peak value as the Lyapunov index a goes from small to large values, although 
a small decrease can be noticed in the a ~ 3A^ region. Finally, from the definition given in 
Eq. ([9]), it is immediate to obtain the average spectral entropy (iij)i which is presented, 
again as function of a, in Fig. [H](c). Its value decreases as the Lyapunov index increases. 
This in turn means that LVs corresponding to smaller positive LEs are more localized in 
Fourier space, i.e. they have more wave-like character than those corresponding to larger 
LEs. However, the decrease in the value of {H^)^ as the Lyapunov index increases is not 
monotonic. We notice that already for a = 200 the average spectral entropy has reached 
its minimum value, but presents a slight increase as a -^ 3A^, which certainly accounts for 
the decrease in the height of the peak in the average spectra of Fig. [6]^a). Our tentative 
conclusion at this point is that the longitudinal LMs are more vaguely defined than in the 
ID case for this density value. 

The results corresponding to the transverse correlation function (C^"-* {k, t)\ are presented 
in Fig. [9l We observe that the a range for which kj^ax has a small value is broader compared 
to the results in Fig [8](a). However, /cmax has a slight increase in value as a ~ 3A^, a 
result which is not consistent with the existence of a transverse Lyapunov mode for this 
a range. This conjecture is supported by the behavior of (C^\kraa,-K-it))^ displayed in Fig. 
[9]^b). This last quantity has its maximum at a ~ 200, with a corresponding minimum of 
(ifj,)^ at the same a value, as seen in Fig. [9](c). This feature of the spectral entropy is 
inconsistent with localization in Fourier space; that is, the wave-like character of the LVs 
is largely diminished in the low a region. Taken together these results make us difficult to 
unambiguously ascertain the existence of transverse LMs in the 3D LJ system. 

For the case of high reduced density p = 0.5 Fig. [TOTa) presents the results for fcmax versus 
the Lyapunov index a corresponding to the longitudinal correlation function (C^^\k^t)\. 
The behavior displayed is very different to that shown in Fig. [8]^a). We first notice that the 
fiuctuations of all quantities are much smaller than those corresponding to the low-density 
value. In the present case the value of fcmax stays close to 1 for a < 100. Then, as a further 
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FIG. 8: (a) Wave number k^ax of the highest peak in the time-averaged longitudinal correlation 
function (^C^\k,t)Y for each a value, (b) The height {C^\kincix-,t)) of the highest peak in 
the time-averaged correlation function, (c) Average spectral entropy {Hj^)^. The parameters are 
p = 0.01, iV = 108, and T = 1.5. 

increases and after a small transient interval, drops rather sharply to a value slightly lower 
than /cjnax = 0.2, which is consistent with that obtained from Fig. [Tl^a). Next, in Fig. 
[TO](b) we observe that the value of {C^\k-[aax,t)). decreases smoothly from its maximum 
at a ^ 3A^ down to its minimum at a ~ 100. Finally, in Fig. [TOlfc) we observe that, at 
variance with the low-density result, the value of the average spectral entropy {Hj^)^ decreases 
monotonically in the whole range of a values, with a minimum for a ~ 3A^ where the LEs 
are the smallest possible ones. From the definition of the spectral entropy we conclude that 
the corresponding spectra for these LVs are most significantly dominated by a few k values. 
Thus the LMs are more sharply defined than in the low density case: a very intriguing state 
of affairs, since no relation was found between the spatial order of the atoms as described 
by S{k) and the highest peak in the LV power spectrum at any density value. 

The results for the transverse correlation function {C^\k,t)\ for the highest density 
studied are presented in Fig. [TTl which display a much reduced fluctuation level, but have 
nevertheless a very similar behavior, compared to those in Fig. [9] corresponding to the low- 
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FIG. 9: Same variables as in Fig. ([8]), but corresponding to the time averaged transversal correlation 
function (C)^\k,t)). Same values of A^, p, and T as in Fig. [H) 




FIG. 10: Same variables as in Fig. ([8]) for the time averaged longitudinal correlation function 
{C^^>{k,t)) corresponding to the high-density value p = 0.5. 
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FIG. 11: Same variables as in Fig. ([8]) for the time averaged transverse correlation function 
{C)^\k^t)\ corresponding to the high-density value p = 0.5. 

density case. That is, the maximum value of {C!^'{k^s,^,t)Y, which is coincident with the 
minimum value of {H^)^, is present in a range of a values that is far from the region in 
which a ~ 3A^, a result not entirely consistent with the existence of transverse LMs. 

V. DISCUSSION 

The overall picture that emerges from our results so far indicates a tangent space dynamics 
much more complicated than that of hard-core systems or the ID LJ fluid. The strong 
fluctuations in the LEs mentioned in Sec. IIIII produce a strong mixture among modes, so 
it is unreasonable to expect that results of the spectral analysis could convey information 
concerning pure modes. Indeed, the maximum wave number fcmax for both the longitudinal 
and transverse correlation functions is highly unstable. In the time interval depicted in Fig. 
[5] for the longitudinal case the /Cmax value very rarely stays close to zero. Rather, it seems to 
wander randomly between zero and one. The instantaneous value of the spectral entropy, 
Eq. ([H]), also seems to have a seemingly random time evolution. Finally, there seems to be 
no correlation between the time evolution of /Cmax and H^°'^ (t) , which was discovered in the 
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case of the ID LJ system and which greatly contributed to a clear-cut definition of the so 
called on and off states [20(]. This is a first indication that the detection of LMs becomes 
more difficult than in the corresponding ID case. 

It turns out that the time averaged longitudinal and transverse correlation functions are 
the relevant variables from which meaningful information about the LMs can be obtained. 
This fact can be understood in terms of the Zwanzig-Mori formalism, in which the Fourier 
components of the fiuctuation of a conserved density vary slowly for a small wave number 



45l | ; from these "slow" variables a meaningful description is then extracted. For the tangent- 
space dynamics the time averaged longitudinal and transverse correlation functions can be 
considered as the "slow" variables. The results presented in Figs. and [7] show that this 
is indeed the case; the spectra are dominated by low wave number values, i.e. fcmax ~ 0.04 
for p = 0.01 and fcmax ~ 0.16 for p = 0.5, both for longitudinal and transverse correlation 
functions, with a corresponding broader peak in the latter case. In both cases /cmax ~ 27r/L, 
which is the smallest nonzero wave number allowed by the periodic boundary conditions used 



(smaller k values are due to the oversampling inherent to the method, see Ref. 46|). Another 
point to be noted is that, for our particular system, the longitudinal correlation function 
reaches its maximum value at a ~ 200, and then its height diminishes slightly, but with the 
peak position fcmax unchanged, as a ~ 3A^. These results can be attributed to a complicate 
mixing of pure modes in the low k regime, which produces the observed degeneracy of the k 
value with respect to the a index. Thus the obtained LMs lose their hydrodynamic character 



and no dispersion relation /cmax vs A'-"-' as those observed in the ID LJ fiuid 20| and hard- 



1 those of 2D 



23|, and even 



core systems J48(] could be detected. Our results even stand in contrast wit. 

coupled map lattices, for which a dispersion relation A ~ fcmax indeed exist 

more sharply with those of the 2D LJ fiuid, for which a corresponding dispersion relation 

has been claimed to hold, although for this last model the evidence seems so far to be 



inconclusive 44 1. 



The most convincing evidence of the existence of longitudinal LMs for the low density 
state comes from the results of Fig. [HI First we observe that the average spectral entropy 
attains its minimum at a ~ 200, a result not entirely inconsistent with those of random 
matrix theory [ij]. Next, as a ^ 3N, the value of this quantity remains close to this 
minimum. Taken together with the already obtained position of the highest peak in the 
power spectrum in the region a ~ 3A^ depicted in Fig. [6](a), these results are compelling 
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evidence of the existence of longitudinal LMs. The fuzziness of the obtained values of the 
reported quantities also led us to suppose that the mixing between modes is strong. Up to 
this point we can affirm that longitudinal LMs do indeed exist, but are more vague than in 
the ID LJ gas and with no hydrodynamical character at all. On the other hand, our results 
on the transverse correlation function displayed in Fig. [HI despite their fluctuations, show a 
tendency that does not allow us to unambiguously classify them as transverse LMs. 

The most important result of this paper was presented in Fig. [10] for p = 0.5. Besides a 
much reduced fluctuation in the studied variables, a very defined jump in k^a.^ at a ~ 100 
is observed. The importance of this fact is that, for the ID LJ gas, a similar jump was 
observed, but for a low-density value [2(1]. In that system this behavior was interpreted in 
terms as of a separation of time scales signaled by a bending in the LS. On the contrary, 
we obtain a sharply defined jump in fcmax at a p value in which the corresponding LS shows 
no bending, as can be observed in Fig. [31 In order to try to understand this seemingly 
puzzling result, we have to remember that, from the suspected importance of hyperbolicity 



14| . the main results for the ID LJ fiuid were obtained 



for the appearance of the LMs 
in a relatively diluted regime 20||, which made this system somewhat similar to hard-core 
systems previously studied 48|. However, for p = 0.5 our system is far from complete 
hyperbolicity, since the combination of attractive and repulsive interactions induces strong 
correlations between the collision events that makes difficult to separate them from each 
other. At the same time, the aforementioned density value is not high enough to make the 
effective interaction among atoms similar to that of a lattice of anharmonic oscillators in 



which hydrodynamic LMs have been detected 25|]. Thus there is no obvious mechanism 
that could account for the sharp jump in k^s.^ when there is no bending in the LS. Therefore 
the role of the alleged separation of time scales, inferred from the LS bending for the ID LJ 
gas, seems to bear no relation to the appearance and phenomenological description of the 
LMs in our 3D LJ fiuid. 

VI. CONCLUSION 

In this paper we have performed a study of the tangent space dynamics of the 3D LJ 
fiuid in order to investigate the possible existence of the LMs which are a distinctive feature 
of the hard-core systems. Our results indicate that longitudinal modes indeed exist for low 
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and high reduced density values, and no conclusive evidence of transverse modes for either 
density studied. The lack of a dispersion relation between the LEs and the maximum wave 
number make both types of modes markedly different from those already encountered in 
other systems. The longitudinal LMs turn out to be much better defined for high values 
of the reduced density. Previously only in hyperbolic systems could LMs be detected at 
high density values [43]. It is highly plausible that by changing the thermodynamic state of 
the system the LMs could present a different behavior than the one reported in the present 
communication. 
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APPENDIX: DYNAMICAL LV CORRELATION FUNCTIONS 

From the Fourier transform of the spatial LV density u^^^ (r, t), Eq. ([S]), we can define 
the intermediate two-time correlation function as 

^(") (k, r) ^ (5(") (k, t + r) S(") (-k, t))t, (A.l) 

which in the 3D case is a second rank tensor with components F^" (k, r). Since the results 
of Sec. IIVCI point clearly to the existence of longitudinal LMs, we will concentrate on the 
longitudinal component of this function. By following the same methodology of Sec. lIVAI we 
are led to an expression for the longitudinal component of the form F^"'^ (kz, t) = (uz {kz, t+ 
t)uz {—kz, t))t. Taking k parallel to the other coordinate directions and averaging we obtain 
the final form F^") (A;, r). For r = we recover the time average of the static LV correlation 
function, i.e. (C*-"^ (/c, t))t = F^' {k, r = 0). The dynamical LV correlation function encodes 
structural as well as temporal correlations, and thus provides more detailed information of 
the system. 

By Fourier transformation with respect to time we obtain the dynamical LV density 
correlation function (DLVDCF) as S*^"^ {k,uj) = (27r)^^ J F^"^ (fc, r) exp(ico'r)(ir. In Fig. 
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FIG. 12: Normalized DLVDCF corresponding to the LV a = 320, with p = 0.01, T = 1.5, and 
N = 108. The curves present results for various values of A: = 27rn/L, n = 1, . . . , 6 from bottom to 
top. The arrow indicates the position uoc of the first peak n = 1. The inset displays the dispersion 
relation u}^"->{k) for different LVs. 



we present the result for p = 0.01 and a = 320 for the longitudinal DLVCF at various integer 
multiples of the wavevector k = 27t/L. It can be observed that, besides the central peak, 
smaller peaks are present as the frequency u increases. This result implies that the tangent 
space dynamics is described by a set of characteristic frequencies, which may have their 
origin in the lack of timescale separation already advanced in Sec. IIIII However, restricting 
our attention to the lowest frequencies, it is clear that the position Uc of the first peak can 
be unambiguously defined for the lowest k value; a similar identification can be made for 
other low k values. Thus a dispersion relation ij^°'\k) for different a values can be extracted. 
The result is presented in the inset of the same figure. It is observed that an approximately 
linear dispersion relation is obtained and thus, besides the characteristic wavevector k^"'"' 
for each LM already determined in Sec. llVCt each mode can also be characterized by a 
frequency ujc{k^°'^). The non- vanishing value of dcu/dk seems to imply propagating wave-like 
excitations. 

The result for the longitudinal DLVDCF at p = 0.5 is presented in Fig. [131 In sharp con- 
trast with the low density state no structure can be observed whatsoever at any wavevector 
number. No characteristic frequency can be easily identified, which implies that no prop- 
agation of any wave-like structures occurs at any wavevector number. Thus the tangent 
space structure described in Fig. (TU] remains unaltered within the time scales studied. More 
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FIG. 13: Normalized DLVDCF corresponding to a = 320, but for p = 0.5. Same T and N values 
as in Fig. [12l 

details of the temporal dynamics in tangent space will be presented in a forthcoming study. 
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